Zero-bias molecular electronics: Exchange-correlation corrections to Landauer's 

formula 
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We show, that standard first principles calculations of transport through single molecules miss 
exchange-correlation corrections to the Landauer formula — the conductance is calculated at the 
Hartree level. Furthermore, the lack of derivative discontinuity in approximations can cause large 
errors for molecules weakly coupled to the electrodes. From Kubo response theory, both the Lan- 
dauer formula and these corrections in the limit of zero bias are derived and calculations presented. 



Much experimental progress has been made in recent 
years in developing methods to measure the conductance 
of single or few molecules in between macroscopic leads, 
and there is keen interest in the theoretical modeling of 
such systems In the case of organic molecules, co- 
valently bound to the metallic electrodes, the transport 
properties are sensitive to the electronic structure, so 
chemical details are likely to be important, and a first 
principles treatment is desirable. In order to describe 
the coupling of the molecule to the macroscopic leads 
in an appropriate way, parts of the leads must be in- 
cluded in the calculation. Given the number of atoms 
required to simulate both the molecule and the first lay- 
ers of the leads, density functional theory (DFT) is an 
obvious choice. 

Since the first sucessful conductance experiments for 
single molecules, there have been several ground-breaking 
calculations of this type, and a variety of codes perform 
DFT-based calculations of I-V curves of single molecules 
between metallic contacts0,0,0,0,0,0,|3. In these cal- 
culations, a potential difference V between the bulk elec- 
trodes is imposed. A self-consistent ground-state Kohn- 
Sham (KS) calculation is performed for the molecule 
plus a few layers of the leads. Then, via Green's func- 
tions the current is calculated using the celebrated two- 
terminal Landauer formula 9]. The macroscopic leads 
enter via self-energies. We denote this "standard" ap- 
proach by its common acronym, NEGF. These calcu- 
lations are parameter-free and often yield qualitative 
agreement with experiment, and so might appear to 
be as rigorous as any DFT calculations. But detailed 
comparison for organic molecules between Au-electrodes 
reveals quantitative discrepancies. Conductances are 
typically overestimated, often by one or two orders of 
magnitude 10]. 

Neither the Hohenberg-Kohn theorem[ll|, which es- 
tablished ground-state DFT, nor the Runge-Gross theo- 
rem for time-dependent nroblems|l2l]. apply to extended 
systems carrying current in homogeneous electric fields. 



In consequence, questions have recently been raised 
about the validity of the NEGF approach [H1E3- For 
example, the calculated transmission is that of the KS 
potential. In the case of a molecule weakly coupled to 
two leads, whose KS levels are sharp, well-separated res- 
onances, the NEGF approach produces peaks in the con- 
ductance at the positions of unoccupied levels of the KS 
system^. Such transitions are known to differ, in gen- 
eral, from the true excitations of the interacting elec- 
tronic system. 

To tackle the transport problem rigorously for a finite 
bias is daunting, and only recently have several sugges- 
tions been put forward |14l Ha . Ha] . In the present paper, 
we examine only the weak bias regime, i.e., the limit in 
which the potential difference across the molecule is in- 
finitesimal, because here we can deduce the exact answer. 

A primary result of this paper is to demonstrate rigor- 
ously that NEGF calculations include only the Hartree 
response of the system. It is alarming, that this level 
of calculation can be inadequate already, when e.g. the 
respnse of an isolated molecule to a static electric field 
is to be considered. This lack is utterly independent of 
which standard approximate functional is used: All DFT 
calculations to date suffer from this limitation. It is in- 
herent in the methodology, just as for all Hartree calcula- 
tions. Second, we estimate the size of the XC corrections 
using the gradient expansion in the current of Vignale- 
Kohn[T3. Even when such contributions are small, the 
lack of derivative discontinuity in semilocal functional ap- 
proximations for the ground state likely produces signif- 
icant errors. Finally, we argue that, under certain con- 
ditions, peak spacings in a zero-bias Coulomb blockade 
experiment are accurately given by NEGF calculations. 

Consider any system that can carry a DC current in 
a specific direction (which we call the z direction) and 
that contains some atomic-sized barrier in this direction. 
For simplicity, we analyze only the symmetric case. We 
apply a weak uniform (also for simplicity) electric field in 
the z-direction, and use time-dependent current density 
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FIG. 1: Double barrier resonant tunnelling cartoon of a 
molecule between two metallic leads; the lowest unoccupied 
molecular orbital (LUMO) has been shifted and broadened 
relative to the isolated molecule. 



functional theory (TDCDFT) to calculate the current re- 
sponse. 

The response of a system to a weak electric field is 
given by the Kubo formula as 



j(rw) = J d 3 r' <7(rr'cj) E cxt (r'cj) 



(1) 



where j(rw) is the first-order (physical) current den- 
sity response to the external electric field, and we use 
atomic units throughout. Here <x(rr'u/) is the frequency- 
dependent non-local conductance describing the response 
to the external field, rather than the response to the total 
electric field. 

Within TDCDFT, the KS system is defined to repro- 
duce the time-dependent current density. Thus Eq. (JTJ 
becomes, for the KS system, 

j(rw) = J d 3 r' CT s [n ](rr'w) (E tot (r'w) + E xc (r'u)) 

(2) 

where j(iw) is still the exact physical current response, 
but now found from the KS non-local conductance (a 
functional of the ground-state density no) applied to the 
KS electric field, which includes both the total electric 
field (external plus Hartree) and XC contributions. 

The external field produces a finite potential drop V 
across the barrier [20j. We restrict ourselves to one dimen- 
sion, to demonstrate the principle. In the limit in which 
lu —>■ 0, the non-local conductance becomes coordinate- 
independent 0, 0| . Thus the left-hand side becomes 
independent of z, and the integral over z' applies only to 
the potentials: 



I = a szz {u; = 0){V + I/ XC ) 



(3) 



where V is the integral over the external and Hartree 
fields and V xc is the induced net XC potential drop in the 
vicinity of the barrier. Equation and the following 
interpretation, are the important results of this Letter. 
We analyze it in two steps. 

(i) Ignore V xc : in the absence of the XC potential drop, 
Eq. 10) tells us that the conductance, I/V, is just that 
of the ground-state KS system. Careful derivations^^ 
show that, for non-interacting systems, 



os Z z{u = 0) = T s (e F )/ir 



(4) 



where T(e) is the transmission through all channels 
through the barrier. The resonances in the KS trans- 



mission function translate into peaks in the conductance 
for the interacting system without correction. 

This brings us to the problem mentioned in the sec- 
ond paragraph, namely the positions of the resonances 
in the NEGF approach compared to the physical system. 
Imagine the case of a one-dimensional double barrier, as 
shown in fig. as the "molecule". Usually, ep is located 
in a spectral gap of the molecule, (as in Fig. QJ, so that 
the system is off-resonance and the conductance is non- 
vanishing only due to the small overlap between the very 
weakly broadened levels and cf- 

To probe the unoccupied resonances at zero bias, ap- 
ply a gate voltage V g to the molecule perpendicular to 
the leads, shifting the LUMO down to ep. As it passes 
through e p (as a function of V g ), there will be a large peak 
in the conductance. But consider what happens when the 
resonance begins to overlap with ep. By virtue of its no- 
dependence, the exact KS ground-state potential differs 
signficantly from the off-resonant case, altering the trans- 
mission characteristics. Peaks in transmission are not at 
the position of (V g = 0)— unoccupied resonances. 

For a sharp resonance, the transmission coefficient is 
given by 



(5) 



where e res and 7 are the position and width of the reso- 
nance which, in DFT, depend on the partial occupation, 
< / < 1, of the resonant level. We will now see how 
the use of smooth, approximate density functionals influ- 
ences the position and width of the resonance. 

Using the spectral function A(E) we can write expres- 
sions for the spectral density of states, n(e) = ^A, as 
well as for the transmission T = %A, to obtain a simple 
linear relationship between n(e) and the transmission of 
such a level, n(e) = 2T(e)/(77r). The self-consistent / is 
found from integrating over n(e) as 



/' = / de n(e) = ~ + -tan" 1 \ 2^ 

2 7T 



.(/) 



1(f) 



(6) 



After inverting this, 

T-\e F ) = 1 + tan 2 {n(f(e F ) - 1/2)} . (7) 

In fig. 2, we plot the transmission over energy for 
this situation, with the parameters given in the cap- 
tion. For this calculation we set the width constant 
(7(/) = 70). The actual dependence of 7 on / is expected 
to be weak and have little effect on transmission peaks. 
NowAe = €h — £l is several eV, where en is the highest 
occupied orbital of the N + 1-electron molecule, and €l 
the LUMO of the iV-electron molecule. As in a NEGF 
calculation using a semilocal functional, e res always de- 
pends smoothly on /, and varies continuously between 
tL and tHi we obtain (assuming e res = ej, + /Ae) the 
dashed line. 
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FIG. 2: Conductance peak from resonance: dashed line is 
self-consistent approximate functional result, dotted line is 
approx result as 7 — > 0, and solid line is exact result. Here 
e L = 0, Ae = 1 and 70 = 0.1. 

For weakly coupled leads, where 7 <C Ae (at any oc- 
cupation), the Fermi level is pinned to the resonance 
(eres(f) -> e F ) for / ^ or 1, so e F = e L + Ae/ and 
we obtain, using Eq. (JJJ the dotted line in fig. 2. Thus, 
in an NEGF calculation, Eq. Q always produces a broad 
peak whose width is comparable to Ae. For the case of 
a linear relation, the width is just Ae/2. 

But this is entirely an artifact of smooth density func- 
tional approximations. The real system has a sharp res- 
onance centered at e#. The exact KS potential of the 
molecule jumps (relative to the reservoir) as soon as there 
is an infinitesimal occupation of the resonant level: 

e r&s = e L + Q{f + r,)&e, (»? — 0). (8) 

Solving eq.Elfor ep, we obtain a peak in the transmission 
of width 7 around e# the solid line in fig. 2. This is the 
famous derivative discontinuity [2l| of DFT. 

Since the true transmission will be much more nar- 
rowly peaked than that in the approximate DFT calcu- 
lation, if the system is off-resonance, the DFT calculation 
produces a strong overestimate of the true conductance. 
Fig. [3 is a cartoon of this situation, in which the width 
of the resonance is 10% of the level shift, and a severe 
overestimate occurs if the Fermi level is at, e.g., 0.5. 

In reality, organic molecules may not be so weakly cou- 
pled to the leads (although the widths of resonances in 
GGA calculations are not a sure indicator of this, for rea- 
sons given above). But this effect, or some remnant of 
it at less weak coupling, would explain the severe over- 
estimate. In Ref. [l(j|, the conductance of benzene was 
calculated in two ways, a) by employing the standard 
approach based on the KS energies and orbitals of an 
equilibrium DFT calculation (GGA) and b) by replac- 
ing the KS data with their counterparts obtained from a 
Hartree-Fock analysis. The typical transmission level in 
a window of 2eV about ep was reduced by an order of 
magnitude. The best current way to see if this effect is 



the culprit for the overestimate would be to perform ex- 
act exchange calculations 22J , which should contain most 
of the effects of the true discontinuity. 

(ii) Include V xc '- Now we discuss how to include: 

V xc = J dz E z , xa {z,u^0) (9) 

where E xc is the XC electric field induced in response to 
the applied field (and ignoring coupling between longitu- 
dinal and transverse modes in a s ). We first note that, 
for any pure density functional, E xc = —V • Sv xc , so 
that Vxc = 3v xc (z — > 00) — Sv xc (z — > —00), i.e., the 
net induced XC potential drop from the extreme left to 
the extreme right of the barrier. In any semilocal ap- 
proximation, Vxc therefore vanishes identically, as the 
induced density response is localized to the region of the 
barrier, so that far from the barrier, Sv xc = 0. Thus, 
using common density functionals, the corrections to the 
KS Landauer formula vanish! 

The corrections are produced by non-local interactions 
present in the exact XC functional. At least two mech- 
anisms are well established which generate such interac- 
tion terms non-local in the density. The most obvious one 
is exact exchange. Even a simple static exchange calcula- 
tion [i^l, including response terms (i.e., beyond NEGF), 
might yield a finite result, i.e., a net drop in the ex- 
change potential across the barrier, just as Hartree does. 
A second mechanism is of the hydrodynamical type and 
therefore finds its natural description within TDCDFT. 
For this case we are able to offer an analytical estimate 
for the size of the effect. In the spirit of Ref. [23J, we use 
the Vignale-Kohn (VK) approximation ^1 to obtain an 
expression for the purely viscous contribution to V xc (al- 
though the original derivation assumes a high-frequency 
regime). It involves a spatial variation of the density, 
n(x), which originates from the backscattering off the 
barrier. A rough estimate is obtained by assuming a) 
that the most important variations in the density profile 
along the wire are of the Friedel-type and b) that the 
viscosity can be approximated by its static, homogenous 
value characteristic of a three-dimensional Fermi liquid. 
With these simplifications eq. I© can be rewritten as 

V^c/V « -(1 - T(e F ))T(e F )/(407r 1 / 2 4 /2 ). (10) 

The viscosity counteracts the current flow and reduces 
the conductance. The factor (1— T) takes into account 
that a barrier causing reflection (and thus density inho- 
mogeneities) is needed for viscous flow to be generated. 

Since kp f» 1, the small prefactor (together with 
the fact that T<1 for well resolved resonances) guar- 
antees only small corrections, as was found in a recent 
calculation j24|. This result, though suggestive, is not 
rigorous, however, since it ignores both the elastic hy- 
drodynamic contribution and the limited validity of VK. 

Finally, we show how, despite the fact that XC cor- 
rections to the voltage, (V xc in Eq. (J3J) are missing, 
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gate charge: Q G = 
Q=4e 
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FIG. 3: Transmission of benzene-dithiolate in the presence 
of a gate. The gate is realized by a square sheet of homogc- 
nously distributed dummy charges (total sum is denoted by 
Qg) located a distance of 3.8 Angstrom above the carbon 
ring. 

NEGF might be used to obtain exact information to be 
compared with experiment. Consider Coulomb blockade 
experiments, which measure zero-bias conductance as a 
function of V g 25]. In Fig. [3] we show the transmission 
of a benzene ring, coupled via two sulphur atoms to gold 
electrodes, obtained from a NEGF transport calculation 
[26j |. The KS LUMO moves towards ep with increasing 
gate charge (voltage). As argued above, the position of 
the LUMO does not give a reliable estimate for the real 
peak position. However, at the particular V g where the 
LUMO passes through ep, its energy must coincide with 
the real many-body level. (In Fig. [31 this would cor- 
respond to Qg ~ 6.) Therefore, at those V g where a 
KS-level crosses ep, a peak in the IV characteristics is 
observed. The peak spacings are given by a ground-state 
DFT calculation. 

Our calculation demonstrates the principle. It is miss- 
ing the derivative discontinuity, but at least we can ig- 
nore the missing V xc - There may be cases where the 
derivative discontinuity is unimportant (i.e., for strong 
coupling or larger molecules) but the missing V xc is not. 
Then standard NEGF calculations will yield accurate re- 
sults for peak spacings in Coulomb blockade experiments, 
although the peak heights are strongly overestimated. 

A not unlikely scenary, in which the Vxc corrections 
are especially large is as follows. Suppose that the qua- 
sistationary non-equilibrium state with flowing current 
can be described by scattering from a single-particle po- 
tential that differs from the ground-state KS potential. 
Then the V xc as appearing in Eq. @ must correct the 
KS off-resonant transmission sufficiently to match that 
of the effective potential. 

We conclude by noting that any formalism to treat a 
many body problem that yields the Kubo response for- 
mula when analyzed within TDCDFT will recover eq. 
(3)0, E|- This work was supported by DOE under 
grant DE-FG02-01ER45928. We thank Roberto Car, 
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